Ab initio quality neural-network potential for sodium 
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An interatomic potential for high-pressure high-temperature (HPHT) crystalline and liquid phases 
of sodium is created using a neural-network (NN) representation of the ab initio potential energy 
surface. It is demonstrated that the NN potential provides an ab initio quality description of multiple 
properties of liquid sodium and bcc, fee, cI16 crystal phases in the P-T region up to 120 GPa and 
1200 K. The unique combination of computational efficiency of the NN potential and its ability to 
reproduce quantitatively experimental properties of sodium in the wide P-T range enables molecular 
dynamics simulations of physicochemical processes in HPHT sodium of unprecedented quality. 

PACS numbers: 31.50.Bc, 31.15.xv, 61.20.Ja, 61.66.Bi 
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I. INTRODUCTION 

Recent experimental studies have shown that sodium, 
a simple metal at ambient conditions, exhibits unex- 
pected complex behavior under high pressure (FIG. H}. 
At ambient conditions, the thermodynamically stable 
form of Na is the highly-symmetric bcc structure, 
which transforms upon compression to the fee phase at 
65±1 GPa 1 and then to a more complex cI16 structure 



at 105±1 GP, 



,2,3 



Above 118 GPa, sodium adopts a 



primitive orthorhombic structure with eight atoms per 
unit cell (oP8), which transforms at 125±2 GPa to an 
incommensurate composite host-guest structure (tI19) 4 . 
It has been reported recently that, at pressure around 
118±2 GPa, sodium crystallizes in a number of very 
complex low-symmetry crystal structures with 50 to 512 
atoms in the unit celli Very recently, a new optically 
transparent double-hexagonal closed packed phase of Na 
(hP4) has been observed above 200 GPa 5 . The be- 
havior of the sodium melting line has also been discov- 
ered to be unusual. Measurements have revealed an un- 
precedented pressure-induced drop in melting tempera- 
ture from 1,000 K at ~30 GPa to room temperature at 
-120 GPaSi 

The complexity of the sodium phase diagram combined 
with experimental difficulties of obtaining detailed char- 
acterization of HPHT phases 6 make molecular dynam- 
ics (MD) simulations an indispensible tool for the inves- 
tigation of the microscopic origins of complex behavior 
in dense sodium. The reliability of MD simulations de- 
pends crucially on the quality of the underlying potential 
energy surface (PES). While density functional theory 
(DFT) provides a comprehensive framework for model- 
ing a wide variety of sodium structures it is not practical 
for lengthy MD simulations because of its high computa- 
tional cost. On the other hand, the construction of accu- 
rate and computationally efficient potentials capable of 
describing various bonding patterns in HPHT sodium is 
a formidable challenge. Many potentials for sodium de- 
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FIG. 1: Experimental phase diagram of sodium 6 . 



veloped based on the pseudopotential theory^ - — and the 
embedded atom model (EAM) 1 ^— are limited to a few 
phases in a narrow P-T region of the phase diagram and 
do not always give a correct description of all properties 
or phenomena of interest. 

In this paper, we present an interatomic potential for 
sodium based on a recently developed high-dimensional 
neural- network (NN) representation of ab initio PESs 18 . 
In this approach, the sodium PES is represented by a 
highly flexible NN optimized to reproduce high-quality 
DFT energies of a large dataset of sodium structures. 
The dataset includes crystal (bcc, fee, cI16, oP8) and 
liquid structures for pressures up to 120 GPa and tem- 
peratures up to 1200 K. We demonstrate that the NN 
potential is capable of reproducing numerous properties 
of sodium phases in this P-T region with an accuracy 
comparable with that of the underlying DFT calcula- 
tions. From a computational standpoint, the NN en- 
ergies, forces, and stress tensor are evaluated with the 
speed of empirical potentials enabling ab initio quality 
MD simulations of sodium on unprecedented length and 
time scales. 
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II. COMPUTATIONAL METHODS 
A. NN representation of ab initio PESs 

Artificial neural networks (NNs), biology-inspired ma- 
chine learning algorithms, are emerging as a new class of 
interatomic potentials that combine the accuracy of an ab 
initio description of PESs with the efficiency of empirical 
force fields^ 8 - - — . In this class, the PES is represented by 
a highly flexible NN capable of describing various bond- 
ing patterns in the system. Given a number of atomic 
configurations for which the ab initio energies are known 
the NN is tuned to reproduce these energies in the train- 
ing process. Overfitting (i.e. obtaining a good fit to the 
training data, but performing less accurately when mak- 
ing predictions) is controlled by testing the performance 
of the NN for an independent test set not used in the 
optimization. Once trained, the NN performs interpola- 
tion to reconstruct the potential energy for new atomic 
configurations with the speed of empirical potentials and 
is, thus, useful to perform lengthy MD runs. 

The NN methodology eliminates many problems asso- 
ciated with empirical potentials. First and foremost, NNs 
completely obviate the problem of guessing a compli- 
cated functional form of the interatomic potential. This 
form is determined automatically by the NN. Second, 
the entire training procedure is fully automated so that 
NNs can be readily extended to new regions of the PES. 
Thus, the significant human effort normally required to 
(re)parametrize the potential is replaced with a short 
computer calculation. Finally, the accurate mapping of 
ab initio energies ensures that all properties determined 
by the topology of the PES are described with the accu- 
racy comparable with that of first-principle calculations. 

Neural networks have been successfully used to in- 
terpolate PESs of simple chemical systems for the last 
decade^—. However, NN-based potentials that can be 
used to map high-dimensional PESs of bulk systems are 
very rar o 18 ' 23 i 27 . Here, we used the NN methodology in- 
troduced recently by Behler and Parrinelloi^. In this NN 
scheme, the total energy of the system is expressed as 
a sum of atomic energy contributions. The atomic en- 
ergies are calculated by a standard feedforward NN as 
a function of the energetically relevant local geometric 
environment of each atom. The local environment of 
a given atom is described by several order parameters 
called symmetry functions, which include radial and an- 
gular many-body terms and depend on the positions of 
all neighbors within a specified cutoff radius. The use 
of symmetry functions (instead of cartesian coordinates) 
as NN input (s) and the partitioning of the total energy 
into atomic contributions ensures that all quantities com- 
puted with the NN (energies, analytical forces, and stress 
tensors) are invariant to translations, rotations, and the 
order of atoms. Futhermore, once the fit is obtained, such 
an NN potential can be applied to systems containing an 
arbitrary number of atoms. 

Details of the high-dimensional NN method are given 



elsewher o 18 ' 23 . Several recent works have demonstrated 
applicability of this methodology to modeling phase di- 
agrams 28 as well as structures of liquids^ and crys- 
tals^. 



B. NN potential for HPHT sodium 

The accuracy of the reference ab initio energies is of 
paramount importance while training the NN. We em- 
ployed the PBE functional in combination with an ultra- 
soft pseudopotential with the 2s and 2p semicore elec- 
trons included explicitly as the valence states. A large 
plane- wave cutoff of 100 Ry and a dense mesh of k- 
points (22 x 22 x 22 for the primitive cell of bcc and 
fee, 12 x 12 x 12 for the primitive cell of cI16 and oP8, 
and 6000k for liquid) were used for all structures so as to 
ensure convergence of the total energy to 1.5 meV/atom. 
The Quantum-Espresso package^ was used to perform 
all ab initio calculations. 

The initial fitting of the sodium NN potential was 
performed on crystal structures that included the zero- 
temperature and randomly distorted structures of bcc, 
fee, cI16, and oP8 phases in the pressure range from -1 
to 200 GPa. Liquid structures of sodium were modeled 
with periodic cubic cells containing 32 and 64 randomly 
arranged atoms with the densities corresponding to the 
0-120 GPa pressure range. The energetically relevant 
local environment for each atom is defined by 48 order 
parameters (see Ref.HHfor details) constructed to include 
all neighbors within 6.0 A cutoff radius. 

After the initial training, the NN was improved self- 
consistently by iterative repetition of the NN-driven MD 
and metadynamics-accelerated Parrinello-Rahman sim- 
ulations 3 ^, collection of new structures emerging from 
the simulations, calculation of the DFT energies for the 
physically relevant structures, and refinement of the NN. 
These iterations were performed until the root mean 
square error (RMSE) of the new structures not included 
in the fit converged to the RMSE of the test set. Af- 
ter the self-consistent procedure, the DFT dataset con- 
tained ^17,000 DFT energies corresponding to more than 
350,000 atomic environments. 10% of all structures were 
randomly chosen for the test set. The best fit was ob- 
tained for a NN with 2 hidden layers, each of which con- 
tains 25 nodes (the total number of the NN parameters is 
1901). The RMSE of the training set is 0.72 meV/atom, 
while the RMSE of the test set is 0.91 meV/atom. 
The maximum absolute errors are 6.18 meV/atom and 
7.62 meV/atom for the training and test sets, respec- 
tively. We would like to emphasize that the fitting proce- 
dure introduces only small error (less than 1 meV/atom) 
in addition to the 1.5 meV/atom numerical (convergence) 
error of the DFT calculations. Thus, the NN-potential 
for sodium is expected to reproduce closely the ab initio 
PES. 
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FIG. 2: Normalized cumulative histogram of the absolute NN 
errors in the training and test sets. It demonstrates that the 
fitting error is less than 1 meV/atom for 83% of structures 
and less than 2 meV/atom for 95% of structures in the test 
set. The inset shows the ordinary histogram for the same 
data. Only a few strucutres in the test and train sets have 
errors larger than 3 meV/atom. 



C. Simulation details 

The zero-temperature structures of sodium crystals 
were obtained by the minimization of the enthalpy with 
respect to the atomic coordinates and lattice parameters 
at constant (external) pressure. The zero-pressure zero- 
temperature bulk moduli and their pressure derivatives 
were calculated by fitting the energy to the Murnaghan 
equation 32 . The second-order elastic constants were de- 
termined from a fit of the energy as a function of an 
appropriate cell distortion to a parabola 33 . 

The lattice dynamics calculations were performed us- 
ing the linear response method within the density func- 
tional perturbation theory 34 . An 8x8x8 g-mesh was 
used in the interpolation of the force constants for the 
phonon dispersion curve calculations. The NN lattice dy- 
namics calulations were carried out with the direct super- 
cell calculation method 3 ^. Convergence tests suggested 
the use of a 6 x 6 x 6 supercell in the force constant cal- 
culation. The PHON program 3 ^ was used to obtain the 
phonon dispersion curves. 

The radial distribution functions g(r), thermal expan- 
sion coefficients, isothermal compressibility, self-diffusion 
and viscosity coefficients for liquid sodium were obtained 
from NN-driven MD simulations. All MD runs were 
performed with the DLPOLY package 3 ^ interfaced with 
the NN code. The temperature was controlled using 
a colored-noise Langevin thermostat that was tuned to 
provide the optimum sampling efficiency over all rele- 
vant vibrational modes 38 . Constant-pressure simulations 
were governed by Nose-Hoover equations of motion with 
Langevin noise on the particle and cell velocitie s 38 ' 39 . 
The time step was set to 1.0 fs. 

State points along several isotherms (T from 400 to 
1000 K with 100 K increment) were obtained from NPT 



simulations with cells of 512 atoms. The density at each 
P and T was obtained by averaging over a 25 ps trajec- 
tory. The volumetric thermal expansion coefficients do 
not change appreciably with temperature and were eval- 
uated for the 800-1000 K range assuming a linear den- 
sity dependence on temperature. The radial distribution 
functions were obtained from 25 ps NVT trajectories for 
systems of 54 atoms. 

To evaluate dynamical properties of a liquid one must 
address the issue of size dependence of the self-diffusion 
coefficient. This effect has been analyzed by Diinweg and 
Kreme r 40 ' 41 who established the following dependence for 
the apparent diffusion coefficient D on the simulation box 
length L: 



D{L) = D(oo) 



2.837297fc B T 



(1) 



where D(oo) is the true diffusion coefficient and r\ is the 
translational shear viscosity, which is much less system 
size dependent than D^. To calculate D(oo) and rj, ap- 
parent diffusion coefficients were computed for different 
system sizes. -D(oo) and r\ were then obtained from the 
y-intcrcept and the slope, respectively, of a linear fit of 
D(L) with respect to 1/L. 

It is important to emphasize that long MD trajectories 
are essential to obtain statistically accurate results for 
transport properties of liquids. Furthermore, it is desir- 
able to perform simulations using large systems. Hence, 
direct ab initio MD simulations (especially with a large 
plane wave cutoff and a dense fc-point mesh) are compu- 
tationally very demanding for the evaluation of the dif- 
fusion and viscosity coefficients 4 ^, whereas, the NN pro- 
vides an affordable and accurate method to determine 
-D(oo) and r\. Each apparent diffusion coefficient was 
calculated from 10 independent 500-ps NVE trajectories. 
Systems containing 256, 512, and 1024 atoms were used 
to determine the dependence of D on L (Figured]). Thus, 
the total simulation time required to obtain D(po) and 
77 at four temperature points is ~60 ns, which clearly 
demonstrates the advantage of the NN approach in com- 
parison with direct ab initio simulations. 



III. RESULTS AND DISCUSSION 

A. Solid phases 

The first test of the NN potential was to calculate 
structural, elastic, and vibrational properties for the 
zero-temperature structures of the two most important 
crystal phases of sodium, bec and fee. The computed 
zero-pressure values for the latice constants and stiffness 
coefficients are summarized and compared with DFT and 
experimental values in TABLE fl] The NN accurately re- 
produces zero-temperature zero-pressure DFT values for 
lattice constants (the error does not exceed 0.02%) and 
all independent elastic constants (the error is less than 
5% for en, C12, C44). The pressure derivatives of the bulk 
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FIG. 3: Diffusion coefficient at T — 490 K as a function 
of the inverse box length. Circles show apparent diffusion 
coefficients with the linear fit represented by the solid line. 
Squares show the diffusion coefficients D(oo) corrected for 
finite-size effects with Eq. [1] The constant dashed line marks 
the value of D(oo). 



TABLE I: Structural and elastic properties of solid sodium 
phases. 







BCC 






FCC 






DFT 


NN 


Exp. 


DFT 


NN 


Exp. 


ao (A) 


4.2008 


4.2018 


4.2908" 


5.2967 


5.2972 


5.4061° 


B (GPa) 


7.63 


7.59 


6.31 a 


7.624 


7.613 


6.433° 


B' 


3.722 


3.501 


3.886 a 


3.71 


3.82 


3.83° 


cii (GPa) 


8.72 


8.70 


8.57 6 


8.74 


8.62 




cu (GPa) 


7.09 


7.03 


7.11 6 


7.07 


7.11 




c 44 (GPa) 


6.11 


6.42 


5.87 b 


5.86 


5.85 




c (GPa) 


0.82 


0.83 


0.73 6 


0.83 


0.75 





a X-ray diffraction study at T = 298 K 1 . 
b Ultrasonic test at T = 80 VM. 



modulus B' are also accurately described. The discrep- 
ancies between the theoretical results and experimental 
data shown in TABLE fl] are due to the fact that ex- 
perimental measurements are taken at non-zero temper- 
atures. When the lattice constants are computed from 
NN-driven MD simulations at 298 K the NN lattice con- 
stant for the bcc structure (4.2816A) is very close to the 
experimentally obtained value of 4.2908A (TABLE Q}. 
The foe lattice constant obtained from MD simulations 
at T = 298 K and P = 70 GPa (3.6570A) is also in very 
good agreement with the experimentally measured value 
of 3.6292A (the foe crystal is not stable at low pressure 
and, thus, we did not attempt to obtain its zero-pressure 
lattice constant). 

FIG. S] displays the calculated pressure dependence of 
selected elastic coefficients. A non-monotonic behavior 
of the tetragonal (d = Cll ~ Cl2 ) and trigonal (C44) mod- 
uli obtained from ab initio calculations is accurately re- 
produced with the NN. Elastic constant softening is a 
usual indication of a pressure-induced phase transforma- 



Pressure (GPa) Pressure (GPa) 
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FIG. 4: Pressure dependence of the elastic coefficients. Points 
(lines) represent NN (DFT) data. 



tion. For example, it has been established that softening 
of the d modulus of bcc is connected with the bcc— ►fee 
structural transition along the tetragonal Bain pat h 44 ' 45 . 
It has also been suggested that the negative melting line 
in the sodium phase diagram could be related to the soft- 
ening of elastic constants in crystal phases^. 

The NN dispersion curves for the bcc and fee struc- 
tures are shown in FIG. [5] for a wide range of pressures. 
Remarkable agreement between the DFT and NN curves 
implies that the NN potential will provide results of DFT 
quality for all finite-temperature properties that are well 
described within the quasiharmonic approximation. The 
NN dispersion curves calculated for bcc at zero pressure 
(FIG.JSJi) are in closer agreement with experimental data 
than those predicted with the embedded atom poten- 
tial o 11 ' 16 and the pseudopotential theory^. We would like 
to point out that the pressure-induced transverse acous- 
tic phonon softening along the [0££] -direction near the 
zone center in bcc is accurately captured with the NN 
(FIG. Eh). This softening is responsible for the instabil- 
ity to the tetragonal deformation (negative d in FIG. HJ) 
and the bcc— ►fee transition mentioned above^ 5 -. 

The calculated enthalpies of sodium crystal phases are 
shown in FIG. [6] Despite small enthalpy diffcrencics be- 
tween the bcc and fee phases (~2 meV/atom) the shape 
of the bcc NN enthalpy curve is in excellent agreement 
with the DFT data. The NN bcc— Kfcc transition pres- 
sure (68 GPa) is only slightly lower than the correspond- 
ing DFT value (77 GPa) and in good agreement with 
experimental measurements (65 GPa at 298 K) 1 . 

Experiments have shown that, in the pressure range 
from 105 to 117 GPa, the cI16 structure becomes the 
stable structure of sodium^— . This structure can be rep- 
resented as a 2x2x2 supercell of bcc with a small shift 
in the positions of the atoms along the space diagonal. 
This shift is described by a single internal parameter x 
(x = gives the bcc structure). Our NN (DFT) calcula- 
tions (FIG. [6] and FIG. [7]) show that cI16 becomes more 
stable than bcc at pressures above 95 GPa (90 GPa). 
The initial rapid growth of the distortion parameter x 
in this pressure region [the dividing line between regions 
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FIG. 5: (color) Pressure dependence of the phonon dispersion curves. Red (black) lines correspond to NN (DFT) frequencies, 
triangles represent experimental data££. 



(i) and (ii)] is accompanied by only a small, on the or- 
der of 0.2 meV/atom, decrease of the enthalpy of cI16 
relative to bcc. Such small energy differencies cannot 
be captured by the NN potential resulting in visible dis- 
crepancies between the x values calculated with the NN 
and DFT. However, this discrepancy is only of minor im- 
portance in MD simulations since fee is the most stable 
phase in this pressure range (i.e. the bcc— >cI16 transition 
is not observed experimentally). It is predicted by both 
NN and DFT that the cI16 strucutre of sodium becomes 
more stable than fee above ~117 GPa at zero temper- 
ature. This transition pressure is only slightly overes- 
timated compared to 105 GPa experimentally observed 
for the fee— »cI16 transition at T = 300 K. Comparison of 
the NN and DFT predictions for the value of the distor- 
tion parameter x demonstrates that the structure of the 
cI16 phase is well reproduced by the NN for the entire 
stability region of cI16 [region (hi) in FIG. [7]. 

The NN and DFT calculations show that the oP8 
structure becomes the stable phase of sodium at pres- 
sures above 160 GPa (the corresponding experimental 
transition pressure is 117 GPa). Thus, the NN poten- 
tial accurately reproduces the experimentally observed 
sequence of stable sodium phases: bcc— >icc— >cI16— S-oP8. 
It predicts the transition pressures in very close agree- 
ment with DFT. The calculated (NN and DFT) zero- 
temperature transition pressures are above the exper- 
imentally observed room-temperature transition pres- 
sures with the discrepancy between theoretical and ex- 
perimental Pt r values increasing with pressure. The pres- 
sure of 150 GPa can currently be taken as the upper limit 
of the validity of the NN potential. Extension of the po- 
tential to higher pressures requires the inclusion of addi- 
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FIG. 6: (color) Stability of crystal phases of sodium relative 
to the fee phase. Solid (dashed) lines represent NN (DFT) 
results. 

tional high-pressure structures into the fitting database. 



B. Liquid 

In this section, we test the performance of the NN 
potential for thermodynamic, structural, and dynamical 
properties of HPHT liquid sodium. 

The temperature and pressure dependence of the den- 
sity of liquid sodium is well described by the NN poten- 
tial. FIG. [8] shows that at zero pressure the NN density is 
in excellent agreement with the experimental and EAM 
data. The pressure dependence of the thermal expansion 
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FIG. 7: Distortion x for cI16 as a function of compression U. 
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FIG. 8: (color online) Thermal expansion of sodium liquid at 
zero pressure. The EAM curve is from Ref. |l7|, the experi- 
mental curve is from Ref. |49 . 

coefficient is shown in FIG. [9] It is well described by the 
Murnaghan equation of state which postulates a linear 
dependence of the bulk modulus on pressure (the line in 
FIG. HD^ 



HID : M„ j 



(2) 



The zero pressure thermal expansion coefficient (an — 
2.861 • 1(T 4 K- 1 ), bulk modulus (B = 5.490 GPa), and 
pressure derivative of the bulk modulus (B' Q — 2.792) 
in FIG. [9] were determined from NN NPT simulations 
at 800 K. It is worth mentioning that the EAM of Be- 
lashchenkc 1 ? predicts negative expansion coefficients at 
high pressures. 

The isothermal compressibility of liquid sodium eval- 
uated from the pressure dependence of the density is 
plotted in FIG. [10] for zero pressure. The NN predic- 
tions are in full agreement with the experimental results. 
The calculated isothermal compressibility at 100 GPa 
(Pt = 3.75 • 10~ 6 MPa -1 ) is an order of magnitude 



FIG. 9: Pressure dependence of the volumetric thermal ex- 
pansion coefficient of sodium liquid. The line is calculated 
using Eq. [5] 
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FIG. 10: Isothermal compressibility of liquid sodium at zero 
pressure. The experimental curve is from Ref. |50| . the dashed 
lines show the uncertanty in the experimental data. 



smaller than the zero pressure value and does not change 
with temperature to any measurable degree. No experi- 
mental data is available for the isothermal compressibility 
at such a high pressure. 

FIG. [TT] shows the radial distribution functions ob- 
tained with the NN potential at two state points (p = 2.8 
and 3.7 g/cm 3 correspond to 50 and 110 GPa, respec- 
tively) . Both of the functions are in good agreement with 
g(r) calculated previously with ab initio MD (p = 2.8 and 
3.7 g/cm 3 correspond to 47 and 107 GPa, respectively)^. 
Ab initio MD predicts slightly overstructured liquid com- 
pared to the NN simulations. The difference between the 
two methods is more noticeable at the higher pressure. 
This discrepancy is most likely due to inaccuracies (low 
plane wave cutoff and/or insufficient fc-point sampling) in 
the ab initio simulations™ The ab initio g(r) in FIG. IllI 
are calculated for a 54-atom system and F-point sam- 
pling of the Brillouin zone (54k points) while the NN is 




FIG. 11: (color online) Radial distribution functions of highly 
compressed liquid sodium. The lower curve is calculated at 
p = 2.8 g/cm 3 and T = 1130 K, the upper curve is obtained 
at p = 3.7 g/cm 3 and T = 930 K and is shifted up by 1. The 
ab initio results are from Ref. [5l] 



trained to reproduce high-quality ab initio energies that 
correspond to 6000k-points. 

The transport properties of sodium liquid at zero pres- 
sure calculated with the NN potential are summarized in 
TABLE [TTJ The experimental values of the self-diffusion 
coefficients are reproduced perfectly with the NN whereas 
the viscosity coefficients are systematically underesti- 
mated. The effect of temperature on the viscousity coef- 
ficients can be described with the exponential law: 



i] = r/ QQ exp(E a /RT), 



(3) 



where E a is the activation energy of an elementary vis- 
cous flow process and rjoa is the pre-exponential factor 
that represents viscosity in the limit of zero activation 
barrier or infinite temperature. Fitting of the data to 
Eq. [3] shows that the activation energies are the same 
for the experimental and NN series (TABLE [TTJ) . It is 
the pre-exponential factor, which includes the activation 
entropy of viscous flow, that is undestimated by ~20% 
in the NN calculations. Nevertheless, the 20% error be- 
tween the calculated and experimental values of r] is con- 
sidered to be very small taking into account that the NN 
potential is derived only from ab initio datas2,. 



IV. SUMMARY AND CONCLUSIONS 

In summary, a NN potential for sodium has been cre- 
ated and tested to reproduce properties of HPHT crystal 
and liquid phases. The NN potential captures the exper- 
imentally observed sequence of pressure-induced solid- 
state phase transformations: bcc^fcc— 5>cI16^oP8. The 



transition pressures as well as structural, elastic, and vi- 
brational properties of bcc, fee, and cI16 phases predicted 
with the NN are in very close agreement with DFT. All 
calculated properties of bcc, fee, and cI16 phases are 

TABLE II: Transport properties of liquid sodium at zero pres- 
sure. 



T (K) 


D (cm 2 ■ sec 


- 1 ) xlO 5 


n (Pa ■ sec 


) xlO 4 




NN 


Exp." 


NN 


Exp. 6 


407 


5.681±0.034 


5.44 


4.709L0.275 


5.810 


448 


7.307±0.033 


7.16 


3.917±0.152 


4.935 


490 


9.070L0.001 


9.05 


3.513L0.003 


4.282 


500 


9.499L0.169 


9.52 


3.388±0.601 


4.152 








NN 


Exp. 




E a (J ■ mol~ L 


) 


6080 


6105 




Voo (Pa ■ sec) 


xlO 4 


0.777 


0.957 



"Ref. 

'Ref. 5^ as reported in Ref. |50 



in quantitative agreement with experimental data. The 
calculated (NN and DFT) zero-temperature transition 
pressures are above the experimentally observed room- 
temperature transition pressures. Further investigation 
of the transition pressure dependence on temperature is 
required to understand the origins of this discrepancy. 
The pressure of 150 GPa (on the NN and DFT scales) can 
currently be taken as the upper limit of the validity of the 
NN potential for crystal phases. The NN potential pro- 
vides an ab initio quality description of thermodynam- 
ics (density dependence on pressure and temperature), 
structural, and dynamical (diffusion and viscosity) prop- 
erties of sodium liquid in the P-T region up to 120 GPa 
and 1200 K. 

The unique combination of accuracy and efficiency of 
the NN potential presented here will lead to dramatic 
enhancement of the quality of MD simulations and bet- 
ter understanding of the microscopic origins of complex 
behavior of HPHT phases of sodium. Detailed studies of 
the nature of the reentrant points on the melting curves 
in the sodium phase diagram and search for new crystal 
phases^, are a few examples of useful follow-on develop- 
ments of this work. 
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